%% 
clear;clc;
%% 

%% input directory of files here
directory = input('Type in full path to directory here, or change in line 6 of this script: ');
cd(    [directory]);
%%

tic
% figure;
for iii = [1:8];

    if iii == 1;
load('KEin_Al2.mat');
    elseif iii == 2;
load('KEin_Al3.mat');
    elseif iii == 3;
load('KEin_Al4.mat');
    elseif iii == 4;
load('KEin_Al5.mat');
    elseif iii == 5;
load('KEin_Ta2.mat');
    elseif iii == 6;
load('KEin_Ta3.mat');
    elseif iii == 7;
load('KEin_Ta4.mat');
    elseif iii == 8;
load('KEin_Ta5.mat');
    end % end iii

disp(['Loaded in KEin  workspaces, ',num2str(iii),' of 8']);

events = length(KE);
%% Determine FOV vs non-FOV particles

% FOV particles will
% 
% 1: Start from z > 13.8713177
%
% 2: Have v velocity (M3) > 0
%
% 3: Be within a circle of radius 3 on the x-y plane when z = 10.97 mm, this is passing through entrance aperture of collimator
%
% 4: Be within a circle of radius 3 on the x-y plane when z = 3.47 mm, this is passing through exit aperture of collimator

%% condition 1
% z_min = 13.8713177;
% 
% vals_FOV_1 = find(z>z_min); 

%% conditions 2

% neg_v_z = find(m3<0);
% 
% vals_FOV_2 = neg_v_z;

%% condition 3
z_enter = 10.97; % mm

xfin_1 = x + m1.*((z_enter-z)./m3); %%% check sign of result (zf-z)
yfin_1 = y + m2.*((z_enter-z)./m3);
zfin_1 = z + m3.*((z_enter-z)./m3);

xyPlane_enter = sqrt(xfin_1.^2 + yfin_1.^2);
vals_xyPlane_enter = find(xyPlane_enter < 3);

vals_FOV_3 = vals_xyPlane_enter;

%% condition 4
z_exit = 3.47; % mm

xfin_2 = x + m1.*((z_exit-z)./m3); %%% check sign of result (zf-z)
yfin_2 = y + m2.*((z_exit-z)./m3);
zfin_2 = z + m3.*((z_exit-z)./m3);

xyPlane_exit = sqrt(xfin_2.^2 + yfin_2.^2);
vals_xyPlane_exit = find(xyPlane_exit < 3);

vals_FOV_4 = vals_xyPlane_exit;

%% determine indicies of particles not in FOV

% a = intersect(vals_FOV_1,vals_FOV_2);
b = intersect(vals_FOV_3,vals_FOV_4);
% c = intersect(a,b);

vals_FOV = b;


allvals = 1:events;
vals_noFOV_1 = allvals(~ismember(allvals,vals_FOV));
vals_noFOV_2 = allvals(~ismember(allvals,vals_FOV_3));

%% Visual diagnostics
% figure;
% histogram(E_D(vals_FOV,1));
% 
% %
% figure;
% histogram(E_D(vals_noFOV,1));
% 
% %%
% dE = 0.025;
% bins = 0:dE:1;
% counts_FOV = histcounts(E_D(vals_FOV,1),bins);
% counts_noFOV = histcounts(E_D(vals_noFOV,1),bins);
% 
% bins_plot = bins(1:end-1) + dE/2;
% 
% figure;
% semilogy(bins_plot(2:end),counts_FOV(2:end),'LineWidth',2);hold on;
% semilogy(bins_plot(2:end),counts_noFOV(2:end),'LineWidth',2);
% 
% r=3;
% teta=-pi:0.01:pi;
% xCircle=r*cos(teta);
% yCircle=r*sin(teta);
% 
% figure;
% plot3(xCircle,yCircle,z_exit-1+ones(1,numel(xCircle)),'b'); hold on;
% 
% for i = 1:100;
% if ismember(i,vals_FOV);
% plot3([x(i),xfin_2(i)],[y(i),yfin_2(i)],[z(i),zfin_2(i)],'Color','g'); hold on;
% elseif ismember(i,vals_noFOV);
% plot3([x(i),xfin_2(i)],[y(i),yfin_2(i)],[z(i),zfin_2(i)],'Color','r'); hold on;
% end
% end
% 
% axis([-10 10 -10 10 -5 20]);
% 
% %%
% A = find(E_D(vals_noFOV,1)>.020);
% Ai = vals_noFOV(A);
%  
% figure;
% plot3(xCircle,yCircle,1+ones(1,numel(xCircle)),'b'); hold on;
% for i = 1:100;
% plot3([x(Ai(i)),xfin_2(Ai(i))],[y(Ai(i)),yfin_2(Ai(i))],[z(Ai(i)),zfin_2(Ai(i))],'Color','r'); hold on;
% end
% axis([-10 10 -10 10 -5 20]);
% 
% figure;
% plot3(xCircle,yCircle,1+ones(1,numel(xCircle)),'b'); hold on;
% for i = 1:100;
% plot3([x(vals_FOV(i)),xfin_2(vals_FOV(i))],[y(vals_FOV(i)),yfin_2(vals_FOV(i))],[z(vals_FOV(i)),zfin_2(vals_FOV(i))],'Color','g'); hold on;
% end
% axis([-10 10 -10 10 -5 20]);
% 
% figure;
% histogram(KE(Ai));

%% free up some memory

clear x* y* z* m1 m2 m3 vals_xy*

%% load in E_D

    if iii == 1;
load('E_D_Al2.mat');
    elseif iii == 2;
load('E_D_Al3.mat');
    elseif iii == 3;
load('E_D_Al4.mat');
    elseif iii == 4;
load('E_D_Al5.mat');
    elseif iii == 5;
load('E_D_Ta2.mat');
    elseif iii == 6;
load('E_D_Ta3.mat');
    elseif iii == 7;
load('E_D_Ta4.mat');
    elseif iii == 8;
load('E_D_Ta5.mat');
    end % end iii

disp(['Loaded in E_D  workspaces, ',num2str(iii),' of 8']);
toc
%% Find % that scatter in and deposit > 20 keV on D1
threshold = 0.020; % MeV
dE = 0.025;
bins = 0:dE:1;

for iB = 1:length(bins)-1;
    A = find(KE>bins(iB));
    B = find(KE<bins(iB+1));
    C = intersect(A,B);
    D_1 = intersect(vals_noFOV_1,C);
    E_1 = find(E_D(D_1,1)>threshold);
vals_noFOV_threshold_1{iB} = D_1(E_1);
    p_scatter_1(iii,iB) = length(E_1)./length(D_1);

    D_2 = intersect(vals_noFOV_2,C);
    E_2 = find(E_D(D_2,1)>threshold);
vals_noFOV_threshold_2{iB} = D_2(E_2);
    p_scatter_2(iii,iB) = length(E_2)./length(D_2);

end
%%
toc
% savename = ['workspace_teeth_20230722.mat']
% 
% save(savename,'p_scatter_1','p_scatter_2');
% disp(['saved p_scatter workspace, ',num2str(iii),' of 8']);

bins_plot = bins(1:end-1)+dE/2;

% plot(bins_plot,p_scatter(iii,:)*100); hold on; 
% 

toc
end % end iii

bins_plot = bins(1:end-1)+dE/2;


save(savename,'bins_plot','p_scatter_1','p_scatter_2');

% 
 figure;
 plot(bins_plot,p_scatter_1(1,:)*100,'b-'); hold on;
 plot(bins_plot,p_scatter_1(2,:)*100,'b--');
 plot(bins_plot,p_scatter_1(3,:)*100,'b:');
 plot(bins_plot,p_scatter_1(4,:)*100,'b.-');
 plot(bins_plot,p_scatter_1(5,:)*100,'r-');
 plot(bins_plot,p_scatter_1(6,:)*100,'r--');
 plot(bins_plot,p_scatter_1(7,:)*100,'r:');
 plot(bins_plot,p_scatter_1(8,:)*100,'r.-');

 
  figure;
 plot(bins_plot,p_scatter_2(1,:)*100,'b-'); hold on;
 plot(bins_plot,p_scatter_2(2,:)*100,'b--');
 plot(bins_plot,p_scatter_2(3,:)*100,'b:');
 plot(bins_plot,p_scatter_2(4,:)*100,'b.-');
 plot(bins_plot,p_scatter_2(5,:)*100,'r-');
 plot(bins_plot,p_scatter_2(6,:)*100,'r--');
 plot(bins_plot,p_scatter_2(7,:)*100,'r:');
 plot(bins_plot,p_scatter_2(8,:)*100,'r.-');